      !<@file Rotinas exigidas por algoritmos de otimização (originalmente o gradiente conjugado do Napack)
      !>@brief Função para avaliação da função objetivo. Retorna somente o valor da função.
      double precision function VALUE(X)

      use globais
!       use globaisCG
      implicit none

      double precision :: X(n) !< vetor de parâmteros da otimização

      integer:: i,j
      integer::iflagF=1 !< indica se será feito o cálculo da função objetivo
      integer::iflagG=0 !< indica se será feito o cálculo do gradiente da função objetivo
      integer::iflagH=0 !< indica se será feita uma aproximação da Hessiana
      double precision :: aux



      write(*,*)'                           _________VALUE()___________'
!       write(*,*)'                           X fcg  '
!       write(*,*)'                         ',X





      call PDMDFDF1(omegas(kw),iflagF,iflagG,iflagH)


      if(kreg.eq.0)then
	  Treg = 0.0d0  !lambdan=0 no BOTH() 
      else if(kreg.eq.1)then
! 	  call regularizaTikho(X,iflagF,iflagG) ! calcula Treg !lambdan  foi lido no arquivo de entrada
      else if(kreg.eq.2)then
! 	  call regularizaL2(X,iflagF,iflagG) ! calcula Treg !lambdan  calculado no BOTH()
      else if(kreg.eq.3)then
! 	  call regularizaWL2(X,iflagF,iflagG) ! calcula Treg !lambdan  calculado no BOTH()
      else if(kreg.eq.4)then
	  call regularizaL2(X,iflagF,iflagG) ! calcula Treg !lambdan  calculado no BOTH()
      else if(kreg.eq.5)then
	  call regularizaL2(X,iflagF,iflagG) ! calcula Treg !lambdan  calculado no BOTH()
      else if(kreg.eq.6)then
	  call regularizaTV(X,iflagF,iflagG) ! calcula Treg !lambdan  calculado no BOTH()
      endif




  
      VALUE = Misfit + lambdan * Treg
      write(*,*) iter,'(iter)',Misfit/value*100.0d0,'%(misfit)',lambdan * Treg/value*100.0d0,'%(reg)'

      write(unitF,100)omegas(kw)/(2.0d0*Pi),iter,Misfit,lambdan,Treg,errorelat*100.0d0
! 
!       C0ant = C0
!       VALUEant = VALUE
    


!       F = 0.0
!       do i=1,n
! 	  F=F+(X(i)-i*10)**2
!       enddo
!       F=F*0.5
!       VALUE = F


      write(*,*)'                           Fobjetivo',value     
      write(*,*)'                           ___________________________'



      return
 100  format(2i10,3E15.6)

      end function

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
      !>@brief Rotina de cálculo do gradiente da função objetivo. Como o cálculo do gradiente depende da simulação (por causa da formulação adjunta), chama-se a rotinha BOTH.
      subroutine GRAD(G,X)

      use globais
!       use globaisCG
      implicit none

      double precision :: X(n)!< Vetor de parâmetros de otimização
      double precision :: G(n)!< Gradiente da função objetivo
      double precision :: VALUE!< Valor da função objetivo

      write(*,*)'                           _________GRAD()____________'

      call BOTH(VALUE,G,X)

      return
      end subroutine



! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
      !>@brief Rotina de avaliação da função objetivo e cálculo de seu gradiente.
      subroutine BOTH(VALUE,G,X)

      use globais
!       use globaisCG
      implicit none

      double precision :: X(n)!< Vetor de parâmetros de otimização
      double precision :: G(n)!< Gradiente da função objetivo
      double precision :: VALUE  !< Valor da função objetivo    


      integer :: iflagF=1
      integer :: iflagG=1
      integer :: iflagH=0

      integer :: i
      integer :: j

      double precision :: xmin,zmin

      double precision :: mapCemX

      double precision :: bestfobjetivo
      double precision :: bestmisfit
      double precision :: bestTreg

      double precision :: Misfitmin
      double precision :: beta =1.0d0
!       double precision :: fobj0
      double precision :: alpha1
      save alpha1
      save beta
      save Misfitmin

      save bestfobjetivo,bestmisfit,bestTreg







      ! solução do problema direto no domínio da frequência
      call PDMDFDF1(omegas(kw),iflagF,iflagG,iflagH)
  



!       if(iter.ne.iterant)then !aproveita p plotar o modelo se for nova iteracao

! 	      xmin = -(Npml-1)*delta
	      xmin = 0.0d0
	      zmin = 0.00d0
! 	      write(*,*)'iplot:',iplot, 'kw:',kw,'iter:',iter, 'Nx0:',Nx0,'Nz0:',Nz0,'N:',N


	      call systempsimage(c,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'veloX',&
		& 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
		& wbox,Lz0/4.d0,Lx0/4.d0, &
		& 0.0d0,Lz0,0.0d0,Lx0,0)!modelo do pd sem PML
! 
! 	      call systempsimage(c,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'veloX',&
! 		& 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
! 		& wbox,Lz0/4.d0,Lx0/4.d0, &
! 		& 0.0d0,Lz0+(Npml-1)*delta,-(Npml-1)*delta,Lx0+(Npml-1)*delta,0)!modelo do pd com PML


! 	      call systempsimage(c0,nx0,nz0,deltaX0,deltaZ0,deltaZ0/2.0d0,deltaX0/2.0d0,'veloX',&
! 		& 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
! 		& wbox,Lz0/4.d0,Lx0/4.d0,&
! 		& 0.0d0,Lz0+(Npml-1)*delta,-(Npml-1)*delta,Lx0+(Npml-1)*delta,0)!modelo do pi, sem a pml


	      ! Salva a aproximacao em cada iteracao para plotagem da evolucao da aproximacao ao longo do processo
! 	      do j=1,Nz
! 		do i=1,Nx
! 		  creal4(i,j) = SNGL(C(i,j)) !retorna single precision malha do PD
! 		enddo
! 	      enddo
! 	      write(unitsol,rec=iplot+1)((creal4(i,j),j=1,Nz),i=1,Nx)!rec=iplot+1 pq iplot começa em zero

	      do j=1,Nz0
		do i=1,Nx0
		  creal4(i,j) = SNGL(C0(i,j)) !retorna single precision malha do PI
		enddo
	      enddo
	      write(unitsol,rec=iplot+1)((creal4(i,j),j=1,Nz0),i=1,Nx0)!rec=iplot+1 pq iplot começa em zero

	      iplot = iplot +1
	      !fim da plotagem da solucao nesta iteracao
!       endif






      if(iterant.eq.-1)then
	  Misfitmin = Misfit
      endif

      if(kreg.eq.0)then !sem regularização
	  lambdan = 0.0d0
	  Treg = 0.0d0  

      else ! com regularização



	  if(iterant.eq.-1)then
! 	      lambdan = lambdanaux*Misfitmin/Treg !primeiro lambdan é calculado aqui, os demais no lambda lbfgs


	      lambdan=misfit

	      treg = 1.0d0
	      do i=1,n
		  gradtreg(i) = 0.0d0
	      enddo


	  else
	      if(kreg.eq.5)then 

		  call regularizaL2(DX,iflagF,iflagG,beta)

	      else if(kreg.eq.6)then 

		  call regularizaTV(DX,iflagF,iflagG,beta)

	      else
		  stop !código de regularização errado
	      endif
	  endif

      endif !





      
      VALUE = Misfit + lambdan * Treg


      do i=1,N
	  G(i) = GradMisfit(i) + lambdan * gradTreg(i)
      enddo

      write(unitF,100)omegas(kw)/(2.0d0*Pi),iter,Misfit,lambdan,Treg,errorelat*100.0d0

      if(iter.ne.iterant)then ! só no inicio das novas iteracoes      
	  write(unitFgrad,100)omegas(kw)/(2.0d0*Pi),iter,Misfit,lambdan,Treg,errorelat*100.0d0
      endif


      iterant = iter

!fim da  Misfit e Grad pela formulação adjunta


! ! ! ! Misfit e grad por diferencas finitas
! ! !       F = VALUE(X) !aproveito só a fobj
! ! !       ! ao final desta iteracao, atualiza valores de F(Xn) e lamda_n
! ! !       FXn = Fobj
! ! !       Xn = X
! ! ! 
! ! !       F = VALUE(X)!aproveito só o Treg
! ! !       lambdan = Fobj/Treg
! ! ! 
! ! !       F = VALUE(X)! agora sim o custo atualizado
! ! ! 
! ! ! 
! ! !       do j=1,N
! ! ! 	  xaux = X(j)
! ! ! 
! ! ! 	  X(j) = xaux + 0.01d0
! ! ! 	  F1 = VALUE(X)
! ! ! 
! ! ! 	  X(j) = xaux - 0.01d0
! ! ! 	  F2 =  VALUE(X)
! ! ! 
! ! ! 	  G(j) = (F1-F2)/0.02d0
! ! ! 	  X(j) = xaux
! ! !       enddo
! ! ! ! fim Fobj e grad por diferencas finitas










!       F = 0.0
!       do i=1,n
! 	  F=F+(X(i)-i*10)**2
!       enddo
!       F=F*0.5
! 
!       do i=1,n
! 	G(i) = (X(i)-i*10)
!       enddo

      if(Misfit.lt.Misfitmin)then
	  Misfitmin = Misfit
      endif

!       if(value.lt.bestfobjetivo)then
! 	  bestfobjetivo =value
! 	  bestTreg = Treg
! 	  bestmisfit = Misfit
!       endif
      write(*,*)'                funcao objetivo:',VALUE,'=',Misfit,'+',lambdan,'*',treg
      write(*,*)lambdan*Treg/Misfit,'(reg/misfit)',lambdanaux,'(peso reg/misfit)'   
!       write(*,*)'                           grad:'   
!       write(*,*)'                         ',G     	
!       write(*,*)'                           ___________________________'



      return
 100  format(1E15.6,1i10,4E15.6)

      end subroutine


! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
      !>@brief Rotina de pré-condicionamento. 
      subroutine PRE(Y,Z)

      use globais
      implicit none

      double precision :: Y(n) !< resultado da multiplicação do vetor Z pela matriz de pré-condicionamento
      double precision :: Z(n) !< vetor a ser multiplicado pela matriz de pré-condicionamento

      integer i,j,kk




!       write(*,*)'                           _________PRE()_____________'
!       write(*,*)'                           OLD',Z
   
!       do i=1,n
! 	  Y(i)=0.0d0 ! zera vetor resultado
!       enddo
!       do i=1,Hnnz
! 	    Y(Hrowind(i))=Y(Hrowind(i)) + Hvalue(i) * Z(Hcolind(i))
!       enddo

!       do i=1,N
! 	  Y(i) = 0.0d0
! 	  do j=1,N
! 	      Y(i) = Y(i) + Hprec(i,j)*Z(j)
! 	  enddo
!       enddo



      do i=1,n
	  Y(i)=Z(i) ! vetor Y = vetor Z  antes de precondicionamento (ou sem precondicionamento)
      enddo

      !Y = H^-1 * Z ==== H Y = Z 




      nrhsH = 1
      ldyH = N


      ! atencao: estrou entrando com Z em Y. Ja sai Y da sol do sistema
      ioptH = 2! aproveitar a matriz ja fatorada para resolver o sistema
      call c_fortran_dgssv( ioptH, N, Hnnz, nrhsH, Hvalue, Hrowind, Hcolptr, &
		  & Y, ldyH, factorsH, infoH ) !CCS  --- no c_fortran_zgssv estou fazendo H^T . Y = Z, mas nao tem problema porque H eh simetrica
      if (infoH .eq. 0) then
	  write (*,*) 'H Solve succeeded'
      else
	  write(*,*) 'H INFO from Solve = ', infoH
	  stop
      endif


	

!       write(*,*)'                           NEW',Y      
!       write(*,*)'                           ___________________________'


      return
      end subroutine


